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ABSTRACT 

We investigate a series of steady-state models of galaxy clusters, in which the hot intracluster gas 
is efficiently heated by active galactic nucleus (AGN) feedback and thermal conduction, and in which 
the mass accretion rates are highly reduced compared to those predicted by the standard cooling flow 
models. We perform a global Lagrangian stability analysis. We show for the first time that the global 
radial instability in cool core clusters can be suppressed by the AGN feedback mechanism, provided 
that the feedback efficiency exceeds a critical lower limit. Furthermore, our analysis naturally shows 
that the clusters can exist in two distinct forms. Globally stable clusters are expected to have either: 
1) cool cores stabilized by both AGN feedback and conduction, or 2) non-cool cores stabilized pri- 
marily by conduction. Intermediate central temperatures typically lead to globally unstable solutions. 
This bimodality is consistent with the recently observed anticorrelation between the flatness of the 
temperature profiles and the AGN activity (Dunn & Fabian 2008) and the observation by Rafferty 
et al. (2008) that the shorter central cooling times tend to correspond to significantly younger AGN 
X-ray cavities. 

Subject headings: conduction - cooling flows - galaxies: clusters: general - galaxies: active - insta- 
bilities - X-rays: galaxies: clusters 



1. INTRODUCTION 

Clusters of galaxies are the largest gravitationally 
bound systems in the universe. They are filled with 
hot gas with T ~ 2 — 10 keV, which loses thermal en- 
ergy prolifically by emitting X-rays. The X-ray surface 
brightness of many galaxy clusters shows a strong cen- 
tral peak that was previously interpreted as the signa- 
ture of a cooling flow with mass accretion ra tes of up 
to several hundred Mq yr^^ fsee lFabianlll994L for a re- 
view). Although the gas temperature is observed to de- 
cline toward cluster centers, recent high-resolution Chan- 
dra and XMM-Newton observations show a remarkable 
lack of emission lines from the gas at temperatures be- 
low about ~ 1/3 of the ambient c luster teinpera ture 
(e-g.- lPetersonet al.ll2001l[2 003: Ta mura et all 120011 : for 
a review see lPeterson fc Fabianii2006l ). In addition, the 
spectroscopically determined mass deposition rates are 
significantly smaller than the classic values estimated 
from the X-ray lumin osity within the cooling regions 
(jVoigt fc Fabianll2004f ). The absence of a cool phase in 
cores of galaxy clusters is suggestive of one or more heat- 
ing mechanisms maintaining the hot gas at keV temper- 
atures for a period at least comparable to the lifetime of 
galaxy clusters. 

Amongst the many candidate heating mechanisms put 
forth recently, there are two leading contenders: 

1. thermal conduction from the hot outer 
regions of the cluster to th e cen- 
ter (e.g., iBertschinger fc MeiksinI Il986t 
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[20mL hereafter ZN03; 



2. heating of the intracluster medium (ICM) by out- 
flows, bubbles, or cosmic rays generated by AGNs 
at cluster centers (e.g., [Briiggcn fc Kaiser 2003, 
Ruszkowski fc Bcgelman 2002, Ruszkowski ct aQ 
200 4 IChaiidran fc Rasera 2007. Guo fc Oh 2008L 



Sijacki et al. 2008). 

Recent theoretical and numerical work (e.g. 
iNaravan fc MedvedevI l2001l : ICho et all [200l) has 
shown that a turbulent magnetic field is not as efficient 
in suppressing thermal conduction as previously thought. 
In particular, iNarayan fc Mcdvedev (2001) showed that 
the effective thermal conductivity k in a turbulent MHD 
medium is a substantial fraction (~ 1/5) of the classical 
Spitzer value Ksp if magnetic turbulence extends over 
at least two decades in scale. On the other hand, 
recent work shows that bouyancy instabilities could 
pote ntially strongly s uppress conductivity in the cluster 
core ()Quataertll200l iParrish fc QuataertI 120081 ) . a point 
we discuss in ^ Following this work, ZN03 shows that 
the electron density and temperature profiles of half of 
the clusters they investigated can be fitted by a pure 
conduction model with the conductivity suppression 
factor / = k/ksp ~ 0.2 — 0.4. However, if only thermal 
conduction operates to balance the cooling, extreme 
fine-tuning of the conduc tion suppression facto r / is 
requ ired ( Guo fc OhI l2008t also see Brcgman fc DavidI 
Il988l) : if / is too low, then a strong cooling flow 
develops, while if / is too high, the temperature proflle 
becomes nearly isothermal, in contrast to observations 
of cool core clusters where the temperature invariably 
declines toward the cluster center. Furthermore, al- 
though thermal conduction is well known to stabilize 
short-wavelengt h pe r turba t ions against therma l in- 
stability (e.g., iFieldl 119651 : iMalagoH et"aLl I1987D . the 
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pure conduction models of the cool core cluste rs are 
thermally unstable against global perturbations (jSokeil 
I2003t iKim fc Nararonl 120031 . hereafter KN03). Using a 
Lagrangian perturbation analysis, KN03 showed that 
the pure conduction model has one globally unstable 
radial mode with the instability growth (e-folding) time 
of ^ 2 — 5 Gyr. Furthermore, if strong perturbations 
are applied (as would be the case in, for instance, a 
cluster mer g er), t he growth times can be even shorter. 
iGuo fc OhI (|2008) showed that if one started from 
arbitrary initial conditions (rather than an equilibrium 
solution), a catastrophic cooling flow quickly develops 
in a conduction-only model with a moderate level of 
conductivity. 

Fortunately, other sources of heating exist. A par- 
ticularly promising candidate is heating by the central 
AGN, for which observ ational evidence has been g row- 
ing in recent years f see iMcNamara fc Nulsei]|2007l for a 
recent review). A majority (~ 71%) of cool cor e clus- 
ters harbor radio sources at their cluster centers (jBurn^ 
[1990). Following the launch of Chandra and XMM- 
Newton, recent high-resolution X-ray observations also 
indicate that these radio sources are interacting with 
their surroundings and often displace the I CM, produc- 
ing X-ray cavities fc.g.. jFabian et ani2000l[Birzan et al.l 
[2004, Forman et al. 200B^ 

By contrast with conduction, heating by the dis- 
sipation of mechanical energy released by central 
AGNs provides a self-regula t ing feedback mechanism 
(e.g., Ciotti fc Ostriker 2001, Ruszkowski fc Bcgelman 
2002, Brighcnti fc Mathews 2003, Kaiser fc Binncy 2003, 
Guo fc OhI [200l . If AGN activity is triggered by 



cooling-induced gas accretion toward cluster centers, 
the AGN heating increases until it halts further accre- 
tion. Thus, the gas accretion rate is self-regulated as 
brief bursts of AG N activity alternate with cooling (e.g., 
IVoit fc Donahuell2005[ l. Due to this AGN feedback heat- 
ing, the accretion flow may automatically adjust itself to 
a low value of the accretion rate, which depends mainly 
on the feedback efficiency e (see equation [T^ . and, in 
a time-averaged sense, the ICM may reach a quasi- 
equilibrium state. This h a s been clearly demonstrated by 
iRuszkowski fc Begelmaiil (|2002l hereafter RB02) in hy- 
drodynamic simulations, who showed that a model clus- 
ter heated by a combination of thermal conduction and 
AGN feedback does not suffer from the cooling catas- 
trophe, but instead r elaxes to a stable quasi-equilibrium 
state. More recently, iGuo fc Otil ()2008D proposed a new 
model of AGN feedback heating, where the ICM is effi- 
ciently heated by both thermal conduction and the cos- 
mic rays produced by accretion-triggered AGN activity. 
In their model, the ICM also relaxes to a stable steady 
state with the mass accretion rate highly reduced, and, 
more importantly, their results do not require fine tuning 
of the various adjustable parameters, including thermal 
conductivity and the AGN heating efficiency. Moreover, 
unlike the conduction-only case, the simulation relaxes to 
a stable state independent of the initial conditions. Al- 
though the detailed dynamics of how the released AGN 
energy is transferred into thermal energy of the ICM may 
be much more complicated than these models and is still 
poorly understood at the present time, these simulations 
strongly indicate that AGN feedback heating may poten- 
tially solve the fine-tuning problem associated with the 



pure conduction model. 

These simulations also suggest that AGN feedback 
heating plays a key role in su ppressing g lobal thermal 
instability in the ICM (see Ro sner fc Tuc ker 1989 for a 
local analysis). While local thermal instability may only 
produce small-scale structures (e.g., local mass dropout, 
emission-line filaments) in galaxy clusters, global ther- 
mal instability may result in a cooling catastrophe and 
a strong cooling flow. Thus, a successful model for the 
ICM must be globally stable, or at least only have in- 
stabilities which grow on extremely long timescales. In 
the present paper, we will use the Lagrangian pertur- 
bation method to formally investigate thermal insta- 
bility in quasi-equilibrium galaxy clusters with thermal 
conduction and AGN feedback heating. Global stabil- 
ity analysis is a method complementary to numerical 
simulations as it allows for the quick identification of 
global trends, quick systematic parameter search and 
helps to build physical intuition. For the spatial dis- 
tribution of AGN feedback heating, we ado pt the ana - 
lytically tractable model proposed bv iBegelmanI (|2001h . 
This model has been corn pared with observations by 
iPiffaretti fc Kaastral (|2006[ ). who show that the model 
usually provides a satisfactory explanation of the ob- 
served structure of cool core clusters, although in a fair 
fraction of their sample the model provide relatively poor 
fits. However, our emphasis is not on the background so- 
lutions of this particular model but their global stability. 
With slight modification, our methods can be applied to 
any model of AGN feedback heating. 

We will show that the feedback mechanism can indeed 
effectively suppress global radial thermal instability in 
cool core (CC) clusters, provided that the AGN feed- 
back efficiency is larger than a lower limit (see § 13.31 and 
§ 13.41 for details). We will also study the dependence 
of the cluster stability on the background ICM profiles 
(§ 13. 5p : for non-cool core (NCC) clusters, which have rel- 
atively flat temperature profiles and which are less stud- 
ied in the literature, the stabilizing effect of the feedback 
mechanism becomes small, but thermal conduction may 
completely suppress global thermal instability. Thus, we 
propose that thermal stability of the ICM favors two dis- 
tinct categories of cluster steady state profiles: CC clus- 
ters stabilized mainly by AGN feedback and NCC clus- 
ters stabilized by thermal conduction. Interestingly, X- 
ray observations suggest that clusters can be subdivided 
into two distinct categories according to the presence or 
absence of a cool core (e.g., Peres et al. 1998, Bauer et al. 
2005, Sanderson et al. 2006, Chen et al. 2007). Our sta- 
bility analysis thus naturally explains these two distinct 
cluster categories. Our model is also consistent with the 
observation by Rafferty et al. (2008) who show that the 
short central cooling time corresonds to younger AGN 
(i.e., shorter X-ray cavity ages) and the anticorrelation 
between the flatness of the temperature profiles and the 
AGN activity recently reported by Dunn fc Fabian 2008. 

The issue of the formation and evolution of NCC and 
CC clusters has recently been addressed by Burns et al. 
(2008) who performed large scale cosmological simula- 
tions. They found out that the CC and NCC clusters 
follow different evolutionary tracks, with CC clusters ac- 
creting more slowly over time and growing enhanced cool 
cores via hierarchical mergers. In contrast, they argued 
that NCC suffered early mergers that disrupted embry- 
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onic cool cores. However, this pioneering work does not 
include the effects of AGN to stop catastrophic cooling in 
the centers of CCs and the numerical resolution in their 
simulations is still too low (15.6/i~^ kpc) to accurately 
study the stability and structure of the cores once they 
are formed. McCarthy et al. (2008) invoked different 
preheating histories to explain the difference between CC 
and NCC clusters. Like us, they find that AGN heating 
is required to stabilize CC clusters. Our present calcu- 
lations make the new suggestion that the bifircation be- 
tween CC and NCC clusters emerges naturally, from the 
fact that clusters with intermediate central temperatures 
are globally unstable. 

The rest of the paper is organized as follows. In § [21 
we describe the time-dependent equations of the ther- 
mal intracluster gas and construct a series of steady- 
state cluster models, in which the mass accretion rate 
is highly suppressed compared to that predicted by stan- 
dard cooling flow models. We then carry out a detailed 
formal linear stability analysis of local and global modes 
of thermal instability in steady-state galaxy clusters in 
§ [3l where we also study the dependence of the clus- 
ter stability on the AGN feedback efficiency and on the 
background cluster profiles. We summarize our main re- 
sults in § |4] with a discussion of the implications. The 
cosmological parameters used throughout this paper are: 
r^TO = 0.3, r^A = 0.7, h = 0.7. We have rescaled ob- 
servational results if the original paper used a different 
cosmology. 

2. STEADY-STATE MODELS 

2.1. Time- dependent equations 

In our model, the intracluster medium is subject to ra- 
diative cooling, AGN feedback heating and thermal con- 
duction. The governing hydrodynamic equations are 
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where k is the effective isotropic conductivity. Depending 
on the details of plasma magnetization and MHD turbu- 
lence, heat transport in the ICM may be very complex, 
and both electron conduction and t urbulent mixing may 
contribute to heat transport (e.g., lLazariaiill2006l ): ad- 
ditionally various instabilities could alt er the nature of 
conductivity within the c ooling region (jQuataertJ 120081 : 
iParrish fc Quataertll2008D . Since at present there is no 
consensus on the nature of conductivity in a turbulent 
magnetized plasma, we adopt the same assumption of 
Spitzer conductivity (with a factor / due to magnetic 
field suppression) that most authors do (e.g., ZN03, 
KN03), 



« = fKSp, 



(5) 



where Ksp is the classical Spitzer conductivity (jSpitzeH 

mil), 
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with the usual Coulomb logarithm InA ~ 37. In this 
paper, we assume that / (0 < / < 1) is constant in both 
space and time. 

According to the ideal gas law, the gas pressure is re- 
lated to the gas temperature T and the electron number 
density n^ via 



pksT A^e , „ 



pm^ 



A* 



(7) 



where ks is Boltzmann's constant, m^ is the atomic mass 
unit, and p and pe are the mean molecular weight per 
thermal particle and per electron, respectively. As in 
ZN03, we use p = 0.62 and p^, = 1.18, corresponding to 
a fully ionized gas with hydrogen fraction X — 0.7 and 
helium fraction Y = 0.28. 

In equation ^ , we neglect the self-gravity of the gas 
and any dynamical effects of magnetic fields. The gravi- 
tational potential <& is determined by the dark matter dis- 
tribution pdm, which we ass ume has a modified N avarro- 
Frenk- White (NFW) form (jNavarro et al.lll997[ ) with a 
softened core (ZN03): 



where d/dt = d/dt -(- v • V is the Lagrangian time 
derivative, p is the gas density, P is the gas pres- 
sure, V is the gas velocity, $ is the gravitational po- 
tential, 7 = 5/3 is the adiabatic index of thermal gas, 
pC = nlk{T) ^ 2.1 X 10"27„2yi/2 gj.gg pj^-3 g-i jg 
the volume coo ling rate due to thermal bremsstrahlung^ 
(jRvbicki fc Lightman 1979. : KN03), and F is the conduc- 
tive heat flux 
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Mo/27r 
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where Tc is a softening radius which introduces a core in 
the dark matter distribution in the very inner regions, Ts 
is the standard scale radius of the NFW profile and Mq 
is a characteristic mass. The corresponding gravitational 
potential is: 



-kVT, 



(4) 



^ We have simplified the form of the cooling function here, ig- 
noring the contribution of metal line cooling, which is important 
at l ower temperatures. Fits t o the full cooling function do ex- 
ist l|Sutherland fc Dopita|[l993I V but they complicate the analytic 
derivation of the global stability analysis. We have experimented 
with the full cooling f unct ion and didn't find qualitative changes 
to our results (see Fig. Illl l. 
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where G is the gravitational constant. The parameters 
Mq and Ts for a given cluster are obtained from the ob- 
served temperature in the outer regions of the cluster. 
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as described in ZN03. We adopt their calculated values 
of these two parameters and their best-fit values of re 
directly. 

The term Ti in equation ([3]) is the volume heating 
rate due to AGN feedback. We ad opt the "efferv escent 
heating" mechanism proposed bv Begelmanl (j2001l ) to de- 
scribe the energy deposition into the ICM by the rising 
bubbles, which are produced by the central AGN. Be- 
cause of the non-negligible gas pressure gradient in the 
ICM, the bubbles will expand as they rise, doing pdV 
work and converting the internal bubble energy to kinetic 
energy of the ICM. The resulting disorganized motion of 
the ICM is quickly converted to heat. Assuming that this 
heating mechanism reaches a quasi-steady state, the de- 
tails of the bubble filling factor, rise rate and geometry 
should can c el. If the cavity expands adiabatically (see 
iGuo fc Otil (|2008f ) for an alternative scenario), the lumi- 
nosity passing through the surface of a sphere at radius 



E cxpb(r) 
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where Pbir) is the pressure of buoyant gas inside bubbles, 
7f, ~ 4/3 is the adiabatic index of buoyant gas (assuming 
it is primarily composed of relativisitic plasma) , and f is 
the unit vector along the radial direction. Assuming that 
the bubble rises subsonically so that pressure equilibrium 
is maintained, Pbii") — Pir), where P{r) is the thermal 
pressure of the ICM, we may rewrite equation (|10p as: 



E 



iagn ( p^ 
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where Pq is the gas pressure at the cluster center, Lagn 
is the AGN mechanical luminosity, and [3 — {-fb ~ l)/7fc- 
AGN activity is likely to be intermittent on a timescale 
of order the Salpeter time ts ^ 10^ yr, and possibly 
as short as ti ^ 10"* — 10^ yr (jRevnolds fc Begelmanl 
[l997). Note that the bubble rise time is typically compa- 
rable to (at most several times) the sound crossing time 
isc ~ lO^'^iooc^ioooy'^ ^'^^ ^ radius r ~ lOOrioo kpc and 
sound speed Cs ^ l OOOcs.ioookms"^ (e.g., see table 3 in 
iBirzan et al.l (|2004f )). and is usually shorter than the gas 
cooling time. Thus, it is justifiable to treat AGN heating 
in a time-averaged sense and assume that the mechanical 
energy of central AGN is injected into the whole ICM in- 
stantaneously (e.g., RB02; Brighcnti & Mathews 2003). 
We further assume the AGN mechanical luminosity to 
be 
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where e is the kinetic efficiency of AGN feedback, and 
Mjn = Airrf^pQVo is the mass accretion rate at the inner 
radius rin, where po and vq are the density and radial 
velocity of thermal gas at Tin, respectively. Therefore, 
the volume AGN heating rate Ti, may be written as 
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where rg is the inner heating cutoff radius, which is de- 
termined by the finite size of the central radio source 
(jRuszk owski fc Begelman"2002': also see a discussion of 
rg in iRoychowdhurv et al. 2005). In the rest of this pa- 
per, ro is taken to be 20 kpc, unless otherwise stated. 

2.2. Steady-state models 

In this subsection, we will construct steady-state clus- 
ter models, which will be used as initial unperturbed 
states in stability analysis presented in the next section. 
We assume that the cluster is spherically symmetric and 
time independent. Equations ([T]), ([2]) and ([3]) are thus 
simplified to 
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where v is the radial gas speed and F is the radial heat 
flux. 

Equations (H]), P^ and pB|) are three ordinary differ- 
ential equations for the three variables P{r), T{r) and 
r'^F{r). We solve these equations as a boundary value 
problem between r — rjn and rout, where we impose the 
boundary conditions: 



ne(nn) = "0, T{rin) = TJ^ 



T{rout) = To 



,F{n, 



(17) 



The first three conditions are taken from either the best- 
fit models of ZN03 or the observational data (see Table 
[T|), while the last condition ensures that there are no 
sources or sinks of heat at the cluster center. Equations 
O, (fT5|) . (I16|) and boundary conditions p7|) form an 
eigenvalue problem with either /, e or M as the eigen- 
value, provided that the other two are given as free model 
parameters. The results of the steady-state cluster mod- 
els presented in this subsection are not sensitive to the 
specific choices for the value of tjh and Tout! here we 
choose Tin = 1 kpc and rout = 1000 kpc, unless otherwise 
stated. We note that the cluster stability does depend on 
the value of rin; this is degenerate with the choice of Tin, 
which is explored in § 13.51 We have experimented with 
different values of rout for our models (e.g., rout is taken 
to be 300 kpc in models of A2597; see Tabled]), and find 
that the results of stability analysis are not sensitive to 
this choice. This is consistent with the extremely long 
gas cooling time in the cluster outer regions. 

To illustrate our results clearly, we need to define 
several physical quantities for each steady-state cluster 
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TABLE 1 

Parameters and results of the steady-state models for typical cool core clusters 





rp a 


q-i a 

J^out 


no" 




M 






r, ll 
'iout 




Name 


(kcV) 


(kcV) 


(cm-3) 


Model 


(Mo/yr) 


e 


f 


(cm 


-') 


h^m/Lx " 


A1795 


2 


7.5 


0.053 <* 


Al 


-0.15 





0.27 


1.42 X 


10-^ 













A2 


-0.15 


0.1 


0.12 


1.26 X 


10-* 


0.52 










A3 


-0.05 


0.3 


0.12 


1.26 X 


10-* 


0.52 


A2199 


1.6 


5 


0.074 


Bl 


-0.015 





0.43 


4.22 X 


10-5 













B2 


-0.015 


0.05 


0.36 


4.00 X 


lO-'' 


0.12 










B3 


-0.00375 


0.2 


0.36 


4.00 X 


lO-'' 


0.12 


A2052 


1.3 


3.5 


0.035 <= 


CI 


-0.006 





0.31 


1.78 X 


10-5 













C2 


-0.006 


0.05 


0.20 


1.62 X 


10-5 


0.29 










C3 


-0.0015 


0.2 


0.20 


1.62 X 


10-5 


0.29 


A2597 ^ 


1 


4 


0.07 9 


Dl 


-0.32 





1.30 


1.39 X 


10-^ 













D2 


-0.32 


0.05 


0.40 


1.31 X 


10-^ 


0.68 










D3 


-0.16 


0.1 


0.40 


1.31 X 


10-^ 


0.68 










D4 


-0.17 


0.1 


0.30 


1.30 X 


10-^ 


0.76 










D5 


-0.14 


0.1 


0.50 


1.33 X 


10-^ 


0.60 



" These boundary values are adopted from the best-fit models of ZN03, unless otherwise stated. 

'"'^out is the corresponding model electron number density at the outer boundary r^Mt- 

"^hagn/Lx is the ratio of the overall AGN heating rate to X-ray luminosity of the cluster. 

'^ Adopted from the Chandra observation I Ettori et al.|[2002). 

■^ Adopted from the Chandra observation ( Blanton et al.ll2001h . 

^To provide a better fit to observations, we take the outer boundary of the cluster A2597 to be rout 

AGN heating cutoff radius to be rg = 40 kpc^ 

8 Adopted from the Chandra observation IIMcNamara et aI1l20011) . 



300 kpc and the 



model. We first define the volume-integrated AGN heat- 
ing rate as 



^agn = / Anr^Tidr, 



(18) 



and the X-ray luminosity as 



Lx = 



Anr pCdr. 



(19) 



A cluster without any heating source will lose its ther- 
mal energy by emitting X-rays. We can thus define the 
isobaric cooling time from equation ([3]) as: 



'^cool — 



7 



1 



p 



(20) 



Table [T] lists the physical parameters and results of the 
steady-state models for four typical cool core clusters. 
Note that the mass accretion rates in our steady-state 
models are usually much less than the Eddington rate 
Med « 26(Mbh/lO^M0)(r?/O.l)-i Mq/jt, where Mbh is 
the mass of the supermassive black hole at the cluster 
center and 77 is the radiative efficie ncy of AGN accre- 
tion. We take the cluster Abell 1795 (jEttori et al.ll2002f l 
as our fiducial cluster. In the first two models (Al and 
A2), the values of e and M are given and / is obtained 
as the eigenvalue. Model Al is a pure conduction model 
(e — 0), while model A2 is a typical hybrid model with 
both thermal conduction and AGN feedback heating. In 
model A3, the AGN feedback efficiency e is chosen to 
be much larger than that in model A2. We further as- 
sume that, in addition to the boundary conditions (J17p . 
the electron density at the outer boundary is fixed (we 



choose the same as that in model A2). Thus, with the 
extra boundary condition, both M and / can be solved 
as eigenvalues. The radial steady-state profiles of Abell 
1795 are presented in Figured] As can clearly be seen, 
the electron density and temperature profiles of these 
three models fit the observational data quite well. The 
profiles of models A2 and A3 are virtually the same, since 
the gas density and temperature at both the inner and 
outer boundaries are exactly the same for these two mod- 
els. Using these cluster profiles and equation ((20| . we can 
then calculate the radial profiles of the isobaric cooling 
time, which are shown in the lower left panel of Figure 
[TJ Obviously, tcooi is much less than the age of the Uni- 
verse within ~ 100 kpc from the cluster center, suggest- 
ing that the radiative cooling is dynamically important 
in the central regions of the cluster. We show the relative 
importance of AGN heating and thermal conduction in 
the lower right panel of Figure[TJ For our hybrid models, 
while thermal conduction is significant in the outer parts 
of the cluster cool core, AGN heating clearly dominates 
at the center. Thus the gas temperature profile in the in- 
nermost regions (< 10 kpc) of the cluster is flatter than 
that in the pure conduction model, as clearly seen in the 
upper right panel of Figure [T] 

Since we will find later that the dependence of the clus- 
ter stability on e varies somewhat with cluster proper- 
ties, particularly the density and tem perature profiles , 
we choose the cluster Abell 2199 (jJohnstone et al.ll2002D 
as our second fiducial cluster and plot its radial steady- 
state profiles in Figure [D The gas cooling time in its 
central cool core (r < 100 kpc) is clearly less than the 
Hubble time. 

ZN03 found that, for some cool core clusters, 
conduction-only models may require implausibly high 
values of thermal conductivity. We apply our model to 
one of these clusters, Abell 2597. As shown in Table 
[Tl the conduction-only model (Dl) for A2597 requires 
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Fig. 1. — Electron number density {upper left), temperature (upper right), isobaric cooling time (lower left), and relative importance 
of AGN heating and conduction [lower right) in three typical steady-state models of the cluster Abell 1795. The dot-dashed line in the 
lower left pane l shows the age of the universe (13.5 Gyr for the cosmology used in this paper). Crosses in the upper panels correspond to 
Chandra data llEttori et al.ll200^ '). See text and Table [l] for additional information. 



a thermal conductivity of / = 1.3. By including AGN 
feedback heating, the required thermal conductivity may 
be much smaller (e.g., / = 0.4 in model D3). Figure [3] 
shows radial profiles of electron temperature and num- 
ber density in our models for A2597 (model Dl - D5), 
which all fit the data reasonably well. We note that 
some cool core clusters may not be well explained by 
mod els with conduction and AG N "effervescent" heating 
(see iPiffaretti fc Kaastral[2006( ). However, they may be 
explained by models with a more realistic AGN heating 
function, and our method of global stability analysis is 
generally applicable to any steady-state AGN feedback 
models with only slight modifications due to the new 
AGN heating function. 

For each cluster, the first model in Table [T] is the pure 
conduction model, where the cooling is entirely balanced 
by thermal conduction. This model determines the max- 
imum value of / for each cluster: as the level of AGN 
heating increases, the required thermal conductivity de- 
creases. In the current paper, we only present our re- 
sults from some typical models with specific values of 
/ (e.g., / = 0.12 in model A3 for the cluster A1795), 
but the conclusions drawn are general to models with 
di fferent levels of thermal co nduction (also see Table 3 
of IPiffaretti &: Kaastrall2006f) . As an example, we addi- 
tionally present models with various values of / for the 



cluster A2597 (models D4 and D5) in Table [T] and Figure 
[3] (also see Figure [7]). 

We note that a minimum level of thermal conduction is 
usually required in our steady-state models, since H/pC 
is not uniform with radius, as readily seen in the lower 
right panels of Figure [1] and [H In other words, the "ef- 
fervescent" AGN heating mechanism (equation [13]) alone 
could not offset cooling at all radii throughout the clus- 
ter. Nonetheless, as we have shown, a combination of a 
moderate level (/ > /min, where /min ^ 0.1 for A1795 
and ^ 0.2 for A2597) of thermal conduction and the "ef- 
fervescent" AGN heating can balance cooling throughout 
the entire cluster, and the resulting steady-state profiles 
of electron density and temperature fit the observational 
data quite well. As / decreases, the heating from thermal 
conduction decreases, and thus the required AGN heat- 
ing increases. When / < /min, the required AGN heating 
in steady-state models with the boundary conditions (|17p 
usually overheats the cluster central regions, while the 
gas temperature and entropy start to drop rapidly with 
increasing radius, which is not consistent with X-ray ob- 
servations. Note that the requirement on conductivity 
(/ ^ /min) may be alleviated if other forms of AGN 
heating (e.g., viscous dissipation of so und waves, see 
Ruszkowski et al.ll2004l : shock heating, seelBriiggen et al.l 
20071 : cosmic ray feedback. see lGuo fc Ohll2008f) are taken 
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into account, since a substantial fraction of sound waves, 
shocks or cosmic rays produced at the cluster center may 
be dissipated in the outer regions of the ICM, with per- 
haps different spatial heat deposition from what we have 
assumed. Further work on this is clearly needed. 

3. RADIAL STABILITY ANALYSIS 

3.1. Perturbation equations 

We linearize equations ([1]) — ([4]) by using the La- 
grangian perturbation method. The background ICM 
is assumed to be in steady state, as described in i; 12.21 A 
Lagrangian perturbation, denoted by an operator A, is 
related to an Eulerian perturbation S in the usual way. 



Here we only consider radial perturbations and 

1 d 



V-^ 



dr 



(r'a 



(24) 



where ^ = Ar denotes the radial component of ^. To 
derive the pertu rbation equations, we use t he following 
properties of A (jShapiro fc Teukolskvl[l983l ): 



A = 5 + ^ • V, 



(21) 



d d 
di~ di ' 



.9 d , dS, d 

A — = — A ^ — . 

dr dr dr dr 



(25) 



(26) 



where ^ i s the dis placement vector of a fluid element 
(see, e.g.. lShapiro" fc Tcukolsky 1983, pp. 130-147). By 
perturbing equations ([T]) and ([7]), we find 



Ap = -pV • t 



AT 
T ^ 



(22) 



(23) 



Applying A to equations ([2]) — (|4]) and using equations 
@ and dUl) - (Uni), we obtain 



dt'^ p dr p dr \ T J p dr dr dr^ 



d_ /Ar\ _ /7 AT _ ae ^ 2e 

dr\T J ^ {2 T dr^ r 



ALr 
47^7-2 



.(28) 
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Fig. 3. — Electron temperature (upper) and number den- 
sity (lower) in five typical steady-state models of the cluster 
Abell 2597. The data points are from the Chandra Observation 
llMcNamara et al.ll2001l) . See Table [T] for additional information. 
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9) 



where Ct = djC/dT\p, Cp = dC/dp\T, U = -Anr^F is 
the radial heat luminosity, and A7i is the perturbation 
of AGN heating rate (equation [T^ 



AH-- 



-r/ro 



ro 1 - e~''/''o 



in- 

r 



-^^AP-H^ 
dP/dr dr dr 



-n 



(/'-D^ 



■f3 



APo 



Po 



AH 



feed; 



(30) 



where APq is the value of AP at the inner boundary Hn, 
and ATifeod = HAM(rin)/A/in is the perturbation of the 
AGN heating rate due to the feedback mechanism, where 
AM(rin) is the perturbation of the mass accretion rate 
at r ^ Tin, 



AM(ri, 



vq dt 



(nn), 



(31) 



which can be easily derived from perturbing the defini- 
tion of the mass accretion rate (M = Airr"^ pv) and using 
equations (P^ and ([M)) . Note that we have neglected 
any time delay between central AGN activity and the 



resulting heating of the ICM, which is justifiable since 
both the AGN duty cycle and the bubble rising time are 
usually much shorter than the gas cooling time (see the 
discussion above equation [13]) . 

Taking ^, AT, AL^ as independent variables, we seek 
solutions of equations ((27|) — (|29p that behave as ~ e'^* 
with time. The term ATifocd in equation (|30p then sim- 
plifies to 



'hi rr 

AHiccd = — ^(nn), 



and equations ^7\ — (|29p may be rewritten as 



(32) 
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AH 
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(35) 



In equations ((33)) — ([35)) and hereinafter, we omit e 
from all perturbation variables. 

Equations (|33p — (|35p form an eigenvalue problem, 
which can be solved numerically with appropriate bound- 
ary conditions to find global eigenmodes and the eigen- 
value a. Before considering the global modes, in the fol- 
lowing subsection, we first study small-scale local modes, 
which may not be capt ured in a global stability analysis 
(see the discussion in iBalbus fc Sokerl 1 1989( 1. The na- 
ture of local thermal stability in a stratified system can 
be subtle, and linked to the convective instability of the 
system (Balbus 1988). 

3.2. Local stability analysis of radial modes 

For local stability analysis, we consider local WKB 
perturbations of the form ^^ ^ik,.r+at ^ Here we neglect 
the feedback mechanism of AGN heating, i.e. ATifecd 
in equation (|30p is taken to be zero. This term simply 
affects the overall normalization of heating without spa- 
tial dependence, and is only important in a global anal- 
ysis. The subsonic background flow is ignored as well, 
so that the unperturbed steady-state ICM is nearly in 
hydrostatic equilibrium. In the local approximation, we 
assume that krr 3> 1 (plane-parallel approximation) and 
that the wavelengths of perturbations are much shorter 
than any spatial scale on which the background quan- 
tities vary (e.g., kr ^ 1/Ap, where Xp = {d\nP/dr)~^ 
is the gas pressure scale height in the ICM). To elimi- 
nate high-frequency sound waves from consideration, we 
also assume that \a\ ^ Cgkr, where Cg — \/ P I p is the 
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isothermal sound speed. Therefore, the perturbed dy- 
namical equation of motion (equation I33p simplifies to 
(KN03): 



ikr^ — ■ 



(36) 



(note that in general, AT/T is complex). 

Similarly, the perturbed energy equation (|35p may be 
rewritten as 



{Pa - p^Cp - pC)ikrS, 



Pa . \ AT 



Taking the leading order terms from the perturbations 
of thermal conduction and AGN heating, we obtain 



A 



1 d 
A'Kr'^ dr 



9 AT 
KTkl — , (38) 



/S.H = 



n 



dP/dr dr 



dP 

A—- = -Hikrt 



(39) 



In the second equality of equation (j39p , the perturbation 
of the radial pressure gradient, A{dP/dr), is evaluated 
by perturbing the Euler equation (equation ^) and as- 
suming that a^ ^ krC^/Xp, which corresponds to slowly 
evolving perturbations^. 
Equations (j36p — ((39|) may be combined to give 



J — IH 7 



J_kT 



(40) 



where 



7 



1 pT^ f dC/T \ 



P 



(41) 



is the growth rate of local isobaric thermal instability in 
the ICM without any heating (e.g.. Field 1965; KN03). 
In the absence of any heating source, equation (|40p re- 
duces to cr = tToo, and we thus immediately recover the 

;eneralized Fi eld criterion for isobaric thermal instability 

BalbuslfToM ) 



dC/T\ 



dT 



<0. 



(42) 



For the ICM with C oc pT^^^, the instability criterion 
given by equation ((42|) is easily met, suggesting that local 
radial perturbations grow exponentially with the growth 
time 

* Provided that \cr\ <^ Cskr, cr^ ^ krC^/Xp is guaranteed as long 
as \a\ < Cs/\p (i.e., the growth time of the perturbation is longer 
than the sound crossing time over a pressure scale height). 



= 0"^ 



= 0.64 



Gyr( 



0.05cm" 



ML] 

2kcV/ 



1/2 



(43) 



Equation (|1(I)) confirms the well known result that 
thermal conduction stabilizes short-wavelength pertur- 
bations (e.g., [Field 1965; Malagoh ct al. 1987; KN03). 
More specifically, for the ICM with C oc pT^^"^, thermal 
conduction alone stabilizes perturbations whose wave- 
lengths (A = 2n/kr) are smaller than the critical wave- 
length 



^Ficld 



= 27r ( 



2kT 
\ipC 

25.6 kpc 



1/2 



0.2 



1/2 



Vo.05cm-3y 



fcsT 
2keV 



3/2 

(44) 



Equation (|40)) also clearly shows that an AGN heating 
term as implemented in our model {TL ex dP/dr) reduces 
the growth rate of local thermal instability, although it 
alone could not suppress local thermal instability com- 
pletely (note that (7 - l)H/(7Pcroo) = 2n/ipC < 2/3, 
since H < pC). We note that local thermal instability of 
the ICM may depend on the actual mechanism by means 
of which the AGN mechanical energy is transferred to the 
thermal ICM; we have assumed complete local dissipa- 
tion of the pdV work done by the expanding bubbles. 
Alternative dissipation mechanisms (e.g., sounds waves 
which damp far away, or the heating effect of dispersed 
cosmic rays) could yield different results. Furthermore, 
in our ID model we have assumed isotropic heating by 
bubbles; in reality the angular variation of bubble cre- 
ation will affect local (and global) thermal stability. 

The local stability analysis is not valid for long- 
wavelength perturbations, and a successful model for the 
ICM must be globally stable. KN03 studied the global 
stability of the pure conduction model, and found that 
it is globally unstable with the typical instability growth 
time of '--^ 2 — 5 Gyr; the growth time can be significantly 
shorter if one applies non-linear perturbations. In the 
next subsection, we will perform a global stability analy- 
sis for our steady-state cluster models presented in § 12. 2i 
and show that the feedback mechanism is essential to 
stabilize global thermal instability in cool core clusters. 

3.3. Global unstable modes 

To analyze global stability of the steady-state models 
of a given cluster, we numerically solve equations (j33p — 
(|35p . which are equivalent to four first-order differential 
equations for the four variables ^/r, AT/T, AL^ and 
V • I = 3^/r -I- rd{£,/r)/dr. We solve these equations 
as an eigenvalue problem, where the eigenvalue is the 
growth rate a. Since we have four variables and one 
eigenvalue, we need to specify five boundary conditions. 
Here we choose the same boundary conditions as KN03. 
The three inner boundary conditions are 



1, 



dr \r 



0, ALr = 0, at r = n^. (45) 
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TABLE 2 

TiMESCALES FOR THE CLUSTERS SHOWN IN TABLE [T] 



^cool.O ^00,0 

Name (Gyr) (Gyr) 



/ Model (Gyr) 



A1795 



A2199 



A2052 



A2597 



0.9 



0.6 



1.1 



0.5 



0.6 



0.4 



0.7 



0.3 



0.28 



0.17 



0.14 



0.07 



Al 


3.8 


A2 


3.3, 43.3 


A3 


stable 


Bl 


2.8 


B2 


4.4 


B3 


16.9 


CI 


6.2 


02 


5.9 


C3 


20.0 


Dl 


2.0 


D2 


3.3 


D3 


38.1 


D4 


27.5 


D5 


20.3 



"icool is the isobaric cooling time at the cluster center. 

^t<xi,o is the growth time of the local isobaric thermal insta- 
bility at the cluster center in absence of any heating source 
(see equation I43J . 

'^igrow is the growth time of the unstable global radial mode. 

Emin is the lower limit of the AGN feedback efficiency, above 

which the ICM is effectively (tgrow > t[j) or completely stable. 

'^There are two unstable global radial modes in this model. 




\ 

Model Bl 

Model B2 
Model B3 




Log r (kpc) 



Fig. 4. — Eigenfunctions of the radial unstable modes for the 
steady-state models of Abell 1795 presented in Fig. [T] plotted 
as a function of radius. The solid lines stand for the lone unstable 
mode for the pure conduction model (model Al in Table[T}. Model 
A2 has two unstable modes: igrow = 3.3 Gyr (dotted lines) and 



Log r (kpc) 

Fig. 5. — Eigenfunctions of the radial unstable modes for the 
steady-state models of Abell 2199 presented in Fig. [2] plotted as 
a function of radius. Each model has one unstable radial mode. 

The first condition is a normalization condition, while the 
second condition guarantees that the solutions are regu- 
lar (due to the presence of a {l/r)d/dr{^/r) term in the 
simplified form of equation I33p . Since rjn is not exactly 
zero, the second condition need not hold strictly. We 
have experimented with different values for the second 
condition (for instance, d/dr{^/r) = ^/r^), and find that 
the results are not sensitive to it. This is consistent with 
the fact that the contribution of AGN feedback (equation 
I32[) in the perturbed energy equation (|35p is independent 
of this condition, although the exact value of the instabil- 
ity growth time for each model will be slightly affected 
due to its contribution to the perturbed cooling term 
{PaV ■ ^) in equation ([35| . The last condition demands 
that the perturbed radial heat luminosity is zero at the 
cluster center, which is equivalent to a zero tempera- 
ture gradient there. The remaining two outer boundary 
conditions are set by the requirement that perturbations 
vanish at the outer boundary, which has cooling times 
much longer than the cluster lifetime: 



(46) 



tgr 

modes. 



: 43.3 Gyr (dashed lines). Model A3 has no unstable radial 



C = 0, AT = 0, at r = j'o 



We use the steady-state models constructed in 
as the background states to calculate the eigenmodes of 
global perturbations. Similar to KN03, we first fix a 
and set AT/T to an arbitrary value at r = rin- We can 
then integrate eqautions ([55]) — ([55]) from r = r;,! to 
r = Tout using a Runge-Kutta method. We use the bi- 
section method to update the inner value of AT/T and 
continue iterating until the first outer boundary condi- 
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tion in equation (j46p is satisfied. Finally, we scan a in 
the range (lO'^Gyr)"! < ct < (lO-^Gyr)"!, and use the 
second condition in equation (j46p as a discriminant for 
solutions. 

We first study cigenmodcs with a real positive a, which 
correspond to globally unstable modes with an instability 
growth time igrow = l/o". The results for the steady-state 
models listed in Table[T]are shown in Table 2, where, for 
comparison, we also list the central values of the isobaric 
cooling time feauation l20p and the growth time too (equa- 
tion [43]) of local isobaric thermal instability in absence of 
any heating source. For our fiducial cluster Abell 1795, 
the pure conduction model Al (/ = 0.27) has one un- 
stable mode with growth time igrow = 3.8 Gyr, which is 
consistent with KN03, who found that the equilibrium 
model of A1795 (/ = 0.2) has one unstable mode with 
t^^ow = 4.1 Gyr. For model A2, where the AGN feed- 
back efficiency is e = 0.1, we found two unstable modes 



with t„ 



3.3 Gyr and 43.3 Gyr respectively. Model 



A3 has virtually the same background ICM profiles as 
model A2 (see § ^T^, but a higher AGN feedback efll- 
ciency (e — 0.3). Our calculation shows that model A3 
has no unstable modes. This is a remarkable result, since 
we, for the first time, show from a linear analysis that 
AGN feedback completely eliminates (radial) global ther- 
mal instability. We note that, without introducing the 
feedback mechanism for the AGN heating (i.e., if ATifocd 
in equation [3D] is taken to be zero), model A3 is still 
unstable and has two unstable modes with igrow = 2.7 
Gyr and 52.8 Gyr respectively. We repeated our calcula- 
tions for many different values of model parameters, and 
found that the models without the feedback mechanism 
are always globally unstable, even when the same models 
with the feedback mechanism included are globally sta- 
ble. Thus, the feedback mechanism of AGN heating is 
the key ingredient for stabilizing global thermal instabil- 
ity in cool core clusters. 

For the cluster Abell 2199, the pure conduction model 
(Bl) is globally unstable with igiow = 2.8 Gyr, which is 
consistent with KN03. With both conduction and AGN 
feedback heating included, models B2 and B3 are still 
globally unstable. However, the growth time of the un- 
stable mode in model B3 (e = 0.2) is tgrow = 16.9Gyr > 
tn, which suggests that thermal instability is dynami- 
cally unimportant and thus is "effectively" suppressed. 
Note that, without introducing the feedback mechanism 
for the AGN heating, the growth time of the unstable 
mode in model B3 is much shorter (igrow = 2.2Gyr). 

In Figure [4] and [5] we plot the eigenfunctions of the 
unstable modes in galaxy clusters A1795 and A2199, re- 
spectively. Clearly, the perturbations have largest am- 
plitude in the cluster central regions (r < 10 kpc) and 
decay rapidly with increasing radius, which suggests that 
perturbations reach nonlinear amplitudes much faster in 
central regions than in outer regions. Such global modes 
have also been found by KN03 for equilibrium ICM mod- 
els with thermal conduction only. 

As shown in Table 2, we obtain similar results for the 
other two cool core clusters listed in Table [TJ It is worth 
mentioning that the pure conduction model for the clus- 
ter Abell 2597 requires an implausibly high value of the 
thermal conductivity (/ = 1.3) and is globally unstable. 
By including AGN feedback heating, model D3 requires 
a smaller conductivity (/ — 0.4) and is effectively stable 



0.1 




0.1 



Fig. 6. — Effect of the AGN feedback efficiency on thermal 
stability in typical cool core clusters. For different models of each 
cluster, the values of / and eM are roughly the same (equal to those 
in the second model listed in Table[T] see the text for details). Top 
panel: scaling of M with e. Bottom panel: the dependence of the 
growth time of unstable modes on e for each cluster. Note that 
the clusters are stable or effectively stable when e is greater than 



a lower limit e„ 



The line for A1795 is double-valued since it 



has two unstable modes; both these modes vanish when e > 0.28. 
Here, we assume e is a constant; if, as suggested by observations, 
e (X A/" (with v ~ 0.3 — 0.6), then e„iin will be smaller (see text). 

(ig,„„ = 38.1Gyr, fore ^0.1). 

3.4. Dependence on the AGN feedback efficiency 

From the results of the previous subsection, one might 
conjecture that the ICM is stable (or effectively stable) 
if the AGN feedback efficiency e is greater than a lower 
limit. To fully explore the dependence of the cluster 
stability on the AGN feedback efficiency, we consider 
the steady-state cluster models subject to the following 
boundary conditions 



T{r out) = Tout, ne(?'out) = "out, 



(47) 



where the value of Uaut for each cluster is chosen to be 
that of the second model for that cluster, as listed in Ta- 
ble [1] As described in § 12. 2i such boundary conditions 
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Fig. 7. — Effect of the AGN feedback efficiency on thermal 
stability in A2597. Each curve is plotted as a function of e with 
a fixed conductivity suppression factor /. Top panel: scaling of 
M with e. Bottom panel: the dependence of the growth time of 
unstable modes on e. For different levels of fixed conductivity, the 
cluster always becomes more stable as the AGN feedback efficiency 
increases. 

ensure that the steady-state ICM profiles are almost ex- 
actly the same for models with different values of e (see 
the lines for models A2 and A3 in Figure [T]) ; it is only 
meaningful to study the dependence of the cluster sta- 
bility on e when the background profiles are virtually 
the same. Since the steady state cluster model only has 
three variables P{r), T{r),r'^F{r) (§[221), the five bound- 
ary conditions (|47p can determine two eigenvalues. Al- 
though our hybrid cluster models formally have three pa- 
rameters (/, e and M), they are essentially determined 
by two parameters: /, which determines the level of ther- 
mal conduction, and eM, which determines the level of 
AGN heating, while the subsonic infiow itself (M) only 
has a negligible effect. Thus, for each cluster, the values 
of / and eM are roughly fixed by the boundary condi- 
tions (|T7)) . and M varies with the free parameter e (see 
Figure [6]) . 

For each steady-state cluster model, we then repeat 
our stability calculations and search for unstable modes. 
The results are shown in Figure [51 where the steady- 
state mass accretion rate and the growth time of unstable 
modes are plotted as a function of the AGN feedback ef- 
ficiency. The value of / for each cluster roughly equals to 
that in the second model of that cluster listed in Table [TJ 



The lower panel of Figure [6] clearly shows that the clus- 
ter is stable or effectively stable (tgrow > tn) when e is 
greater than a lower limit emin- We note that this result 
generalizes as well to models with different values of /, as 
seen in Figurc[71 which shows the effect of AGN feedback 
efficiency on global stability for models of A2597 with 
/ = 0.3, 0.4 and 0.5. Assuming that the real intracluster 
gas is in a stable quasi-steady state, our global stability 
analysis thus suggests a constraint on the kinetic effi- 
ciency of AGN feedback. As listed in Table 2, the values 
of Emin for these four typical cool core clusters are emin ~ 
0.07 — 0.28, which is roughly consistent wit h the recent 
estimate of e ~ 0.3 for radio-loud AGNs bv lHeinz et al.l 
()2007( ) and is marginally consis tent with o b servat ional 
estimates of e ~ 0.01 - 0.1 by lAllen et all (|2006f ) and 
iMerloni fc Heind (|2007f) . In our AGN feedback model, 
we assume that e is a constant. However, recent X-ray 
observations seeni to sugges t that e oc M", where i/ ~ 0.3 
(jAllen et al.ll2006D or 0.6 (M erloni fc Heina[2007f) . In this 
case, the AGN feedback will be even stronger (i.e., in 
equation [301 AT^feed = HAM(ri„)/Min -h HAe/e), and 
thus the value of emin (evaluated at M appropriate for 
the steady-state case) will be reduced. For the cluster 
A1795, Cmin is reduced from 0.28 to 0.22 {v = 0.3) or 
0.18 {v = 0.6). For the cluster A2199, emin is reduced 
from 0.17 to 0.13 (i^ = 0.3) or 0.10 {v = 0.6). The ac- 
tual value of Cmin is also sensitive to the exact form of 
the AGN heating law adopted, particularly its spatial de- 
pendence. Here we have only considered the pdV work 
due by rising b ubbles, and igno red, for instance, cos- 
mic ray heating ("Guo fc Olill2008i), vis cous dissipation of 
sound waves JR uszkows ki et al.l |2004() or shock heating 
(jBriiggen et al.l 12007 1. However, while we do not place 
great store in the absolute value we obtain for emin, the 
fact that there is a minimal heating efficiency Cmin for a 
given temperature profile should be fairly robust, since 
it only requires that Lagn oc eM. 

3.5. Dependence on the background profiles 

In this subsection, we will study the dependence of 
the cluster stability on the background steady-state ICM 
profiles. We adopt the clusters AbeU 1795 and AbeU 2199 
as our fiducial clusters. Since the gas cooling time at the 
outer boundary is much longer than the Hubble time, 
we consider steady-state cluster models with fixed val- 
ues of Tout and riout , which are chosen to be the same as 
those in the second model listed in Table [1] for each clus- 
ter. Of the three inner boundary conditions in equation 
(fT7)) . we choose two (including rf^F{rin) = and a given 
value of Tin for a specific model). We have four bound- 
ary conditions for three first-order ordinary differential 
equations (§ 12. 2p . Thus, the value of M can be solved 
as the eigenvalue for a specific steady-state model with 
a given value of / and e (which are deemed to be fixed 
by the physics of thermal conduction and black hole ac- 
cretion respectively) . The corresponding central electron 
number density for each steady-state model (represented 
by the varying value of Tin) is shown in the upper panels 
of Figure [51 and [TPl In addition, it is possible to vary the 
inner boundary rjn. Varying rjn at fixed Tm is basically 
denegerate with varying Tjn at fixed rin, and thus we do 
not explore this additional dependence. 

We first consider the models of the cluster A1795 with 
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Fig. 8. — Dependence of the cluster stability on the background 
profile for the cluster A1795. For a fixed value of / and e, the cor- 
responding central electron number density {upper panel), steady- 
state mass accretion rate (middle panel), and the growth time of 
unstable modes [lower panel) are plotted as a function of the cen- 
tral gas temperature Tin. The dot short- dashed line in the lower 
panel shows the growth time of the unstable mode in steady-state 
models with / = 0.12, e = 0.1 and no feedback mechanism for 
AGN heating (i.e., ATif^od in equation 1301 is taken to be zero), 
while the dot long- dashed line stands for the central gas cooling 
time in models with / = 0.12, e = 0.1. 

/ = 0.12 and e = 0.1 (i.e., variations of model A2). The 
dependence of the cluster stability on the background 
profile (represented by Tin) is plotted in Figure [8] (soHd 
line), which clearly shows that the model with either a 
relatively flat (Tin > 4.5 keV) or steep (Tjn < 1.7 kcV) 
temperature profile is (effectively) stable, while ther- 
mal instability can develop at intermediate temperatures. 
Figure O shows radial profiles of electron number density, 
temperature and the eigenfunction (^/r) of the radial un- 
stable mode for three typical steady state models with 
Tin = 6 keV, 3 keV, and 1 keV respectively. 

3.5.1. Cool core versus non-cool core clusters 

X-ray observations also suggest that clusters can be 
subdivided into two distinct categories according to the 
presence or absence of a cool core (e.g., Peres et al. 1998, 
Bauer et al. 2005, Sanderson et al. 2006, Chen et al. 
2007). In this section we discuss how our model may 
account for this effect. 

The dot-dashed line in the lower panel of Figure [5] 
shows the growth time of the unstable mode in steady- 
state models with / = 0.12, e = 0.1 and no feedback 
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Fig. 9. — Radial profiles of electron number density {upper 
panel), temperature {middle panel) and the eigenfunction of the 
radial unstable mode {lower panel) in three typical steady-state 
models of the cluster Abell 1795 with / = 0.12 and e = 0.1. 

mechanism for AGN heating (i.e., ATifocd in equation [501 
is taken to be zero). Without the feedback mechanism, 
the instability growth time in cool core clusters with rela- 
tively steep temperature profiles is very short ('^ 2 Gyr), 
suggesting that the feedback mechanism plays a key role 
in stabilizing thermal instability in these clusters. As 
the central gas temperature increases and the central 
gas density decreases, the stabilizing effect of the feed- 
back mechanism becomes smaller, which is reasonable 
since the perturbation of the central mass accretion rate 
(note that Tagn oc A/jn) scales as the central gas density 
(equation [311) and the importance of the feedback mech- 
anism in the energy perturbation equation p5p may be 
expressed as ATifocd/ (-Per V • $), which (using equations 
(dSl [HI [311) and dP/dr - pg) scales as T"^ at any spe- 
cific radius. Colder, denser (i.e. low entropy) gas has a 
steeper pressure gradient in a fixed potential well, which 
increases the volumetric rate at which cavities perform 
pdV work as they rise. 

On the other hand, although the stabilizing effect of 
the feedback mechanism becomes negligible for non-cool 
core clusters with relatively fiat temperature profiles, our 
stability analysis surprisingly shows that these cluster 
models are also (effectively) stable. As shown in ij3.21 
the diffusive AGN heating {H ex dP/dr) only increases 
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Fig. 10. — Dependence of the cluster stability on the back- 
ground profile for the cluster A2199. For a fixed value of / and e, 
the corresponding central electron number density (upper panel), 
steady-state mass accretion rate (middle panel), and the growth 
time of unstable modes (lower panel) are plotted as a function of 
the central gas temperature Tin. The dotted line in the lower panel 
stands for the Hubble time, while the dot short-dashed line shows 
the growth time of the unstable mode in steady-state models with 
/ = 0.36, e = 0.2 and no feedback mechanism for AGN heating, 
while the dot long- dashed line stands for the central gas cooling 
time in models with / = 0.36, e = 0.2. 

the growth time of local thermal instability, while ther- 
mal conduction may completely suppress local perturba- 
tions with wavelength less than Apicid, which increases 
as the central gas temperature increases and which may 
then be greater than the cooling radius. Thus, we ex- 
pect that thermal conduction alone may completely sup- 
press thermal instability in these non-CC clusters. This 
is confirmed by our stability analysis: we constructed 
NCC models with only conduction and found them to 
be stable. 

The short and long dashed lines in Figure [S] show the 
models with a higher thermal conductivity and AGN 
feedback efficiency respectively. Obviously, a higher level 
of thermal conduction increases the stability of the clus- 
ter. However, a higher level of AGN feedback efficiency 
only increases the stability of CC cluster models, but 
has a negligible stabilizing effect on non-CC cluster mod- 
els. Thus, AGN feedback heating is not required in NCC 
clusters. On the other hand, for CC clusters with lower 
central temperatures, AGN heating is generally required: 
otherwise one is either unable to build an equilibrium 
profile with a physically plausible level of thermal con- 




FlG. 11. — The growth time of unstable modes for models of the 
cluster A2597 with / = 0.4 and e = 0.05, plotted as a function 
of the central gas temperature. Both the stable CC and NCC 
branches are seen in the calculations with either our simplified 
cooling function or a full cooling function (see text). 

duction / < 1, or else the CC cluster is globally unstable 
on fairly short timescales. We did the same calculations 
for the cluster A2052 and A2597, and found similar re- 
sults. 

We stress that, except for models with very high val- 
ues of / and e (which may always be stable, as seen 
from the trend of different lines shown in Figure [8] and 
[TU|) . the clusters are usually stable either when the cen- 
tral temperature is high (non-CC) or when the central 
temperature is sufficiently low. Figure [8] (and the anal- 
ogous Figure [10] for A2199) suggest that the interme- 
diate values correspond to globally unstable solutions. 
This is consis tent with ID hydrodynamic simulations of 
RB02 and iGuo fc OhI (J2008): if conductivity is large, 
the clus ter relaxes to sta ble NCC states (see Figure 1 
and 2 of lGuo k Ohll2008f) : if / is small (i.e., conduction 
couldn't offset the cooling), the cluster first cools grad- 
ually through NCC states, and then quickly relaxes to 
st able CC states (s ee Fig. 1 in RB02 or Fig. 4 and 5 
of IGuo fc Ohll2008[ ). We note that the radiative cooling 
times in non-CC clusters are generaly longer than in the 
CC ones and in some of the non-CC clusters no heating 
is required to prevent thermal instability on a timescale 
shorter than the Hubble time (note however that in the 
cases shown in Figures |H] and [10] (dot long-dashed lines) , 
the central cooling timescale is always shorter then the 
instability growth time). Irrespectively of whether this 
is the case or not, our model naturally explains the di- 
chotomy between CC clusters where AGN feedback sig- 
natures are observed (and where the AGN play the key 
role in establishing global stability) and the non-CC clus- 
ters, where conduction may play a role but where the 
AGN feedback is not readily observed. This trend for 
the non-CC (flat temperature) clusters to have no AGN 
as opposed to CC clusters with clear AGN observational 
signatures is clearly seen in the recent data compiled by 
Dunn & Fabian (2008). This is also consistent with the 
observation by Rafferty et al. (2008) who show that the 
short central cooling time corresponds to younger AGN 
(i.e., shorter X-ray cavity ages). 

Metal line cooling may become important when the 
gas temperature is low. We checked our calculations 
with a full cooling fun ction (equation 35 of 'Guo fc Obi 
l2008l which is based on [Sutherland fc Dop ita 1993), and 
did not find qualitative changes to our results. As an 
example. Figure [TT] shows the dependence of the cluster 
stability on the central gas temperature for the cluster 
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A2597, which has the lowest Tin in our cluster sample. 
With both free-free and metal line cooling included, the 
gas cooling rate increases and the dependence of igrow 
on Tin becomes similar to that of the higher-temperature 
cluster A1795 (see Fig. l]), but both the stable CC and 
NCC branches still exist. 

Through hydrodynamic numerical simulations, RB02 
shows that the ICM heated by a combination of AGN 
feedback and thermal conduction usually relaxes to a sta- 
ble quasi-steady state. Surprisingly, our stability analysis 
shows that a specific steady-state model may be globally 
unstable if the AGN feedback efficiency is lower than a 
critical value. This "inconsistency" may be explained if 
the cluster with lower e relaxes to a steady state with 
lower central gas temperatures, which may then be effec- 
tively stable, as clearly sh own in this subsec tion. Recent 
numerical simulations bv lGuo fc OhI (|2008D indeed con- 
firms that the cluster central regions in their cosmic-ray 
feedback models with lower e cool to higher densities and 
lower temperatures in the final steady state (see Figure 
7 of IGuo fc Oh (2008)). 

3.6. Global decaying modes 

We have also searched for global decaying modes with 
real and negative a. For each steady-state model of the 
ICM, similar to KN03 (see Figure 4b of KN03), we found 
a series of decaying modes, within which there exists a 
slowest decaying mode with the smallest decay rate. Fig- 
ure[T2l shows the eigenfunction {£,/r) of the slowest decay- 
ing mode in typical steady-state models for the cluster 
Abell 1795 and AbeU 2199. The decay rate of the most 
slowly decaying mode could potentially serve as an indi- 
cator of the "attractor" solution toward which a cluster 
evolves after it has been reset by a merger, and might 
be worthy of further study. It may also be interesting 
to explicitly study the effect of AGN feedback on global 
non-radial modes and overstable modes, which we leave 
to fut ure work (see lMalagoli et aLlll987LlBalbus fc Soked 
Il989l and KN03 for local analyses of these modes in the 
ICM with thermal conduction). 

4. SUMMARY AND DISCUSSION 

Recent Chandra and XMM-Newton X-ray observations 
suggest that the hot intracluster gas is maintained by one 
or more heating sources at keV temperatures for a pe- 
riod at least comparable to the lifetime of galaxy clusters 
(§[!]). Since the emission lin es expected from cooler g as 
are notoriously absent (e.g.. lPeterson fc FabianI (|2006( l). 
the heating source (or sources) should also effectively 
suppress the thermal instability of the ICM. Although 
thermal conduction stabilizes thermal instability for per- 
turbations with short wavelengths, the equilibrium ICM 
heated by thermal conduction alone is globally unstable 
and will evolve to a cooling catastrophe under global per- 
turbations. On the other hand, using numerical hydro- 
dynamic simulations, RB02 and Guo fc Oh (2008) show 
that the ICM, which is initially in a state far from equi- 
librium and is heated by a combination of thermal con- 
duction and AGN feedback, usually relaxes to a quasi- 
equilibrium steady state. Although these simulations 
adopt simplified ID models for the elusive spatial distri- 
bution of AGN heating, they strongly suggest that the 
AGN feedback mechanism plays a key role in suppressing 
global thermal instability. 
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Fig. 12. — Eigenfunctions of the slowest decaying modes (ct < 0) 
in typical steady-state models for the cluster: (a) Abell 1795 and 
(b) Abell 2199. For each steady-state model, there exist decaying 
modes with larger decay rates, which are not shown in this figure. 

In this paper, we perform a detailed formal analysis of 
thermal instability in the ICM with both AGN feedback 
heating and thermal conduction. To build initial un- 
perturbed states for stability analysis, we first construct 
steady-state cluster models, where the gas density and 
temperature profiles fit observations quite well and where 
the mass accretion rates are highly suppressed compared 
to those predicted by the standard cooling flow mod- 
els fi) 12. 2p . Using the Lagrangian perturbation method, 
we then derive a set of differential equations that form 
an eigenvalue problem for global radial modes with the 
mode growth rate a as the eigenvalue. For pure con- 
duction models of typical cool-core clusters, we find that 
the ICM has one unstable radial mode with the typical 
growth time ~ 2 — 6 Gyr (Table 2), whic h is consistent 
with the results of lKim fc Naravcm (|2003l ) . However, for 
the hybrid models with both AGN heating and thermal 
conduction, global thermal instability is effectively re- 
duced or even completely suppressed if the feedback effi- 
ciency e is greater than a lower limit emjn. Interestingly, 
if the AGN heating (equation [T3|) is independent of the 
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central mass accretion rate (i.e., no feedback mechanism 
for AGN heating), the ICM is stiU unstable (§[331), which 
suggests that the feedback mechanism is essential to sup- 
press thermal instability. 

Assuming that the real intracluster gas is in a sta- 
ble quasi-steady state, our global stability analysis thus 
suggests a minimum value emin of the kinetic efficiency 
of AGN feedback, if it is to suppress a cooling flow. 
The value of emin required in typical cool-core clusters is 
around ^ 0.07 — 0.28 (see Table 2), which is roughly con- 
sistent with the estimate of the j et production effic iency 
(~ 30%) for radio-loud AGNs bv lHeinz et all (l200l and 
which is marginally consistent w ith recent observat ional 
estimations of e ^ 0.01 - 0.1 bv lAllen et all (|2006[ ) and 
iMerlo ni fc Heinj ()2007[ ). Note that the value of emin wiU 
be reduced if e is an increasing function of the central 
mass accretion rate as suggested by recent observations 
(§ 13. 4p . Although the existence of emin should be fairly 
robust, its exact value will also depend on the assumed 
form of AGN heating law. 

A related important issue in AGN feedback models is 
how and what fraction of the cooling gas at a distance of 
order 1 kpc gets to the cluster center and finally fuels the 
AGN. This could be a very complex process due to the 
large range of distance scales involved. Our calculation 
extended from ~ IMpc in the cluster outer regions to 
~ Ikpc at our innermost integration point. However, the 
black hole gravitational radius of influence tbh is much 
further in: 



rBH = 5 — = 0.05 kpc f 



Af, 



BH 



VIO^Mq/ VsOOkms 

(48) 
It is conceivable that if the flow is steady all the way 
down to the black hole, it transitions from a cool- 
ing flow to an adiabatic Bondi flow (see, for instance, 
Quataert fc Narayan 20(10 and the discussion in § 5 of 



Chandran fc Rasera ,200^ . and there is tentative ob- 



servational evidence that the Bondi formula may pro- 
vide a reasonable approximati on of th e accretion process 
in X-ray luminous galaxies (jAllen e t al. 2006). How- 
ever, there are a host of potential complications: an- 
gular momentum, which would cause a torus or accre- 
tion disk to form instead, magnetic pressure and outflows 
(jProga fc Begelmanll2003l ). and local thermal instability 
resulting in the formati on of stars and cold gas blobs 
(jPizzolato fc Sokeiil2005D , resulting in "cold accretion" or 
the possibility that the black hole is fed by stellar winds. 
In this paper, in line w ith most numerical simulations 
in the literature (e.g., iRus zkowski fc Begelman' '2002!; 
Brighenti & Mathcw j [2003 : Hoeft & Briiggcn 2004; 
Vernaleo fc Revnoldsll2006[ ) we have assumed that essen- 



tially all of the inflowing gas in our innermost bound- 
ary point eventually makes it to the black hole, over 
timescales long compared to the AGN duty cycle, but 
shorter than or comparable to the gas inflow timescale, 
iflow ~ lO^yr. This behaviour is likely to be intermit- 
tent rather than steady: gas accumulates in an accretion 
disk/torus around the black hole, which eventually trig- 
gers an outburst and leads to gas consumption, etc. We 
assume that in steady state, the black hole will even- 
tually consume all the gas which is supplied to it. In 
particular, in our model, while the black hole can exert 
thermal feedback on the cooling gas (through its heat- 



ing activity), it exerts negligible hydrodynamic feedback 
(by consuming gas faster or more slowly than the supply 
rate, thus affecting pressure forces throughout the cool- 
ing flow region). If indeed such a scenario applies, the 
results of our present paper should be recalculated in de- 
tail, since in that case M cannot be freely varied, but is 
further constrained by the accretion law. Given the large 
uncertainties and difficulty of this calculation, we leave 
this to future work. 

In § 13.51 we study the dependence of the cluster stabil- 
ity on the background steady-state proflles and find that 
the unstable cluster models with e < Cmin may become 
effectively stable if the central gas temperatu re drops to 
much lower values. Numerical simulations bv lGuo fc OhI 
(|2008f l indeed confirm that the cluster central regions in 
their cosmic-ray feedback models with lower e usually 
cool to higher densities and lower temperatures in the fi- 
nal steady state. Thus, unlike the pure conduction mod- 
els, where nonlinear evolution of global unstable modes 
lead to the cooling catastrophe, the ICM in our hybrid 
cluster models with AGN feedback included may always 
evolve to a quasi-steady state, which is effectively stable. 
On the other hand, we also show that thermal conduction 
will completely suppress thermal instability in non-CC 
clusters with relatively flat temperature proflles. Thus, 
the stability of the ICM favors two distinct categories 
of cluster steady state proflles: CC clusters stabihzed 
mainly by AGN feedback and non-CC clusters stabilized 
by thermal conduction. Interestingly, recent X-ray ob- 
servations also suggest that clusters can be subdivided 
into two distinct categories according to the presence or 
absence of a cool core (e.g., Peres et al. 1998, Bauer et 
al. 2005, Sanderson et al. 2006, Chen et al. 2007). 

It is perhaps even possible that these two categories 
of clusters represent different stages of the same object. 
The importance of thermal conduction on global scales 
obviously depends on the large scale structure of the 
cluster magnetic fields. Recent calculations suggest that 
thermal conduction of heat into the cluster core can be 
self-limiting: in cases where the temperature decreases 
in the direction of gravity, a buoyancy stability sets in 
which re-orients a radial magnetic field to be largely 
tra nsverse, sh utting off conduction to the cluster cen- 
ter (jOuataertl [20081 . Non-linear simulations indicate the 
heat flux could be reduce to ~ 1% of the Spitzer value 
(Parrish & Quataert 2003). Thus, the following scenario 
could arise: as conductivity falls, gas cooling and mass 
inflow will increase, triggering AGN activity. The rising 
buoyant bubbles or the gas convective overturn mediated 
by cosmic rays (Chandran 2004) may re-orient the mag- 
netic field to be largely radial again, increasing thermal 
conduction and reducing mass inflow, shutting off the 
AGN until the heat flux driven buoyancy instability sets 
in once again. The cluster could therefore continuously 
cycle between cool-core (AGN heating dominated) and 
non cool-core (conduction dominated) states. 

Global linear stability analysis can clearly serve as a 
useful complement to simulations, both for the physi- 
cal insight they can deliver and speed in exploring pa- 
rameter space. Formally, a globally stable equilibrium 
state is a necessary but insufficient condition for a suc- 
cessful cluster model. A linear analysis fails for large 
non-linear perturbations (as a cluster would undergo, for 
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instance, during mergers). Nonetheless, the linear sta- 
bility analysis qualitatively reproduces the same features 
we have observed in ID hydrodynamic simulations when 
we start the simulation from arbitrary initial conditions: 
conduction-only models suffer drastic cooling flows, while 
conduction -|- AGN heating models relax to a stable state 
with low mass inflow rates. For a given value of e and 
/, a wide range of stable profiles are accessible. From a 
suite of linear stability analyses alone, it is not possible to 
predict the detailed final temperature and density profile 
the cluster relaxes to (which is somewhat sensitive to ini- 
tial conditions), but it is possible to make general state- 
ments about stability and whether the cluster will relax 
to a cool core or non cool-core configuration. Finally, 
we have restricted our attention to a ID radial analy- 
sis. While a 2D or 3D analysis would be more general, it 
would be much more complicated and unlikely to prove 



enlightening, at least for global modes. Non-radial modes 
may prove important in the case of local thermal insta- 
bility, but the fastest growing unstable global mode — the 
possible cooling flow we wish to stem — is likely to respect 
spherical symmetry, unless the heating sources (such as 
the AGN jets) are severely anisotropic. Such models are 
beyond the scope of the present paper. 
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